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ABSTRACT 


-*.  * 


Highly  efficient  algorithms  have  been  derived  for  calculating  the 
probability  of  detection,  P^,  for  several  target  fluctuation  medeis.  These 
include  the  Marcum  model  (constant  target  cross  section)  and  the  four 
Swerling  models  (chi-squared  fluctuations).  Programs  are  presented  in 
computer  free  notation  as  well  as  FORTRAN  (for  accuracies  of  10  ^ and 
10  ^)  for  each  of  the  models. 

In  addition,  an  approximate  expression  for  the  elevation  of  Pp, 
applicable  when  the  target  return  is  constant  during  a scan  but  exhibits  log- 
normal fluctuation  from  scan-to-scan,  is  presented,  together  with  a list- 
ing of  a FORTRAN  program  that  evaluates  the  expression. 

An  efficient  iterative  routine  for  the  determination  of  the  threshold 
from  the  probability  of  false  alarm  is  also  included. 


in 


LIST  OF  ILLUSTRATIONS 


Fig*  Page 

1 (k,  m)  pairs  contributing  to  ^P^fX,  Y). 8 

2 Program  for  evaluating  PL(X,  Yk  11 

3 Approximate  limits  on  X and  Y vscfor  single  precision  ...  17 

4 Program  for  the  Marcum  model.  Case  0..... 18 

5 Program  for  Swerling,  Case  1 ..............  21 

6 Program  for  Swerling,  Case  3 23 

7 Program  for  Swerling,  Case  4.. 26 

8 Program  for  subroutine  PGEN{N,X,Y,  ZK,  P,  YMO) 38 

9 Probability  of  detection  vs  X for  N =’  1,  PpA  = *0  ^ 

(Cases  0-5) .. 40 

10  Probability  of  detection  vs  X for  N = 10,  PpA  = 30  ^ 

(Cases  0 - 3 and  5) . 41 

1 1 Probability  of  detection  vs  X for  N = iOO,  P_A=1Q  ^ 

(Cases  0-3  and  5)  .42 


TABLES 

Page 

A -I  45 

B-I  60 


B-II 


63 


SECTION  1 


INTRODUCTION 


\m 


Some  highly  efficient  algorithms  for  the  evaluation  of  probability  of 

detection  on  N incoherently  integrated  returns  have  been  determined  and  are 

presented  here.  The  efficiency  is,  of  course,  inversely  related  to  accuracy, 

consequently  two  FORTRAN  versions  are  given:  one  with  accuracy  of  10 

-12 

and  the  other  with  an  axcuracy  of  1 0 

In  April  I960,  the  IRE  Transactions  on  Information  Theory  published 
papers  by  Marcum  and  Swerling  [Refs.  1, 2,  3]  that  derived  expressions  for 
the  probability  of  detection,  P^,  for  targets  of  constant  backs cattering 
(Marcum)  and  of  four  different  models  for  backscatter  fluctuation  (Swerling). 

The  resulting  formulas  were  manipulated  into  more  convenient  forms  [Refs,  4, 
5,  6,  ],  and  graphs  supplementing  those  of  Marcum  and  Swerling  were  furnished. 
Heidbreder  and  Mitchell  [Ref.  7]  added  the  log-normal  model  for  characterizing 
the  fluctuation.  This  latter  model  has  been  found  to  be  very  useful  for  targets 
that  have  large  mean  to  median  ratios.  Other  models  exist  [Refs.  8,9, 10], 
but  we  will  limit  ourselves  to  these  since  they  cover  most  situations  well 
enough. 

Although  the  graphs  provided  in  the  references  are  useful,  an  efficient 
means  of  computation  is  sometimes  necessary.  Two  problems  arise  if  one 
uses  the  expression  as  presented.  First,  computations  can  be  quite  lengthy. 


Ill to rial 


especially  if  great  accuracy  is  required,  and  second,  there  can  be  problems 
of  overflow  and  underflow  even  though  all  answers  are  between  zero  and  one. 
New  efficient  expressions  have  beer  obtained  for  the  Marcum  and  Swerling 
models,  and  these  are  modified  to  reduce  the  computation  effort  further 
without  sacrificing  accuracy.  An  approximation  is  presented  for  the  case  in 
which  the  radar  cross  section  is  constant  within  a scan  but  exhibits  log-normal 
fluctuations  from  scan  to  scan. 

Programs  for  the  evaluation  of  the  various  P^fs)  are  presented  in 

computer  free  language  and  in  two  FORTRAN  versions  for  an  IBM  3?0;  the 

•6  m 1 2 
first  with  an  accuracy  of  10”  and  the  second  with  an  accuracy  of  10”  . 


SECTION  2 


EXPRESSIONS  FOR  PROBABILITY  OF  DETECTION 


In  this  section  the  expressions  for  P^  for  the  six  cases  are  presented 
without  proof  since  they  are  derived  in  the  literature  [Refs.  1 , 2,  3,4, 5,  6, 7]. 
In  order  to  preserve  the  Swerling  model  numbers  1 through  4,  we  will  refer 
to  the  constant  target  as  case  0,  the  Swerling  cases  as  cases  1 through  4,  and 
the  log-normal  fluctuation  as  case  5. 

2.  1 Case  0: 

The  Marcum  derived  expression  for  P^  is 


0Pn(X,Y)  = 


(2.1) 


where  N is  the  number  of  pulses  incoherently  integrated 
Y is  the  threshold  level 

X is  the  average  signal-to-noise  ratio  of  a single  pulse 
Xjv  is  the  total  signal-to-noise  ratio  of  all  N pulses 
I (Z)  is  the  nth  order  modified  Bessel  function. 


if  the  N pulses  are  not  identical,  we  still  have 


X = XN/  N . 


By  substituting  in  Eq.  (2. 1)  the  infinite  series  for  I^fZ) 


00 


k=0 


(Z/2)2k 
k!(n  + k)l 


(2.2) 


and  interchanging  the  order  of  summation  and  integration,  one  obtains  Fehlrnr  *’ s 
equation  [Ref.  4] 


0Pn(X.Y)  = c 


y ~ »k  N-l  ~Hc 

■ "Ett  £ «■« 

k=0  m=0 


2.  2 Case  1 : 

The  signal-to-noise  ratio,  X,  with  mean  X,  is  assumed  constant 
throughout  the  integration  scan  but  distributed  as  chi  squared  with  2 degrees 
of  freedom  from  scan  to  scan.  The  resulting  is  [Refs.  5,  6] 


2.  3 Case  2 


The  signal-to-noise  ratio-  X,  with  mean  X,  it,  assumed  to  be  distributed 
as  chi-equared  with  2 degrees  of  freedom  and  independent  from  pulse  to  pulse. 
The  expression  for  PD  is  [Refs.  5, 6] 


2 N (X,  Y)  = 1 


E 

m=N 


Y / Y \m 
1 + X \1  +X/ 


mf 


(2.5) 


2.4  Case  3 

The  signal-to-noise  ratio,  X,  with  mean  X,  is  assumed  to  be  dis- 
tributed as  chi-squared  with  4 degrees  of  freedom  from  scan  to  scan  but 
constant  within  a scan.  The  expression  for  is  [Refs.  5,  6] 

Y _ 

7 / f~f/7\  V \ 

for  N = 1 


1+X/2A  + Jg/21.Y.\ 
\ (1  + X/2)2j 


3^*N  ^ ~ 


N - 1 -Y 
,Y‘  e 


(M  - 2)  I (1  + xn/2) 


N-2 

E 

m=0 


-Y  Ym 

e 1 for  N > 2 

m ! 


+ e 1 +XnA  A +^i-\N'2  [i  - 2Ll±  4 —1— 

\ V2/  L V2  1+V2 


N-2 

i-E e 

m=0 


* + V2 


(nfci) 


m 


m! 


(2.  6) 


2.  5 Case  4: 

The  signal -to -noise  ratio,  X,  with  mean  X,  is  assumed  to  be  distributed 
as  chi-squared  with  4 degrees  of  freedom  and  independent  from  pulse  to  pulse. 
The  expression  for  is  [Refs.  5,6] 


N _ k " _ 

-,/vvi-l  / 1 \V  N!  WV  . 1 + X/2  (1  + x/2)m 


(2.7) 

2.  6 Case  5: 

The  signal-to-noise  ratio,  X,  with  mean  X,  is  assumed  to  be  log- 
normally  distributed  [Ref.  7]  from  scan  to  scan  but  constant  within  a scan. 
The  log-normal  distribution  with  parameters  X and  p is 


(x  I P ) 


•JZv  <r: 


2o-2 


0 S x < 


(2.8) 


A I 

where  M = — is  the  median  of  X and  tr  = v 2i  np  is  the  variance  of  InX.  The 
P 

expression  for  is 


Jf*  oo 

f 0PN(x,Y)px(x|X,p)dx 
0 


(2.9) 
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SECTION  3 


NEW  EXPRESSIONS  FOR  PROBABILITY  OF  DETECTION 


In  this  section  the  expressions  for  cases  0 through  4 are  modified  so 


as  to  achieve  more  efficiently  a specified  accuracy  and  also  to  eliminate  the 


problem  of  overflow  and  underflow.  The  efficiency  is  a result  of  the  3 facts: 


(1)  a recursive  method  is  used  to  evaluate  the  probability,  (2)  the  bound  on 


the  error  consists  of  two  factors,  one  in  terms  of  and  one  in  terms  of  Y, 


each  of  which  is  approaching  zero,  and  (3)  the  terms  used  for  the  bound  are 


similar  to  those  used  in  the  evaluation  and  require  little  additional  calculation. 


3.  1 Case  0: 


We  begin  with  Eq.  (2.  3)  and  change  the  order  of  summation  (see  Fig.  1) 


to  obtain 


QPN(X,Y)  = e 


■(XN  +Y) 


N-l  °°  SF  °°  __  . 00 

Lym  xn  , y y1 

rn  • ' -J  k!  Zj  m *•  * 


m=N 


k=m+l-N 


■Y  Y' 


« / m-N 

m=N  \ k=0 


<*•» 


since 


co 

~ = eX 


(3.2) 


We  separate  Eq.  (3.  1)  into  two  terms,  PL(X,  Y),  where  with  LiN 


N-l 


pl'x’y>=E 


-Y 


m=0 


m=N 


( 


and 


OO 

Rl(X,Y)=  ^ e~Y 
m=L+l 


( 


so  that,  if  we  use  PL(X,  Y)  to  represent  0PN(X,  Y),  our  truncation  error 
will  be  Rl(X,  Y).  We  can  bound  RL(X,  Y)  by 


/ 


L 


/ 


L+l-N  t?  k\ 


If  we  choose  L large  enough  such  that 


(r  L V / y L+l-N  y k\ 

'••',tfr  (-■■  2 irV 

m=0  f \ k=0  / 


(3.6) 


where  € is  our  desired  accuracy,  then  we  have  guaranteed  that 


0Pn(X,Y)  - Pl(X,Y)£€. 


(3.7) 


A program  in  computer -free  notation  to  evaluate  P^(X,  Y)  is  shown 

wAp 

in  Fig.  2.  The  notation  used  is  based  on  one  devised  by  Iverson  [Ref.  11], 


‘The  arrow  notation  within  the  box  symbolizes  a specification,  that  is,  the 
statement  SUM  *-  YMS  is  translated  to  mean  chat  the  quantity  SUM  is  specified 
by  YMS.  A branch  is  denoted  by  an  arrow  outside  the  box  leading  to  the  next 
statement  to  be  executed.  A comparison  is  denoted  by  a colon  (:),  and  a 
branch  i3  executed  if  the  comparison  condition  specified  on  the  arrow  is 
satisfied;  otherwise  the  next  instruction  in  sequence  is  executed.  As  an 
example,  we  interpret  the  segment  below  as  follows:  if  M = N-l,  go  to  line 
n+5,  otherwise  set  M to  M+l.  Multiply  YM  by  Y,  divide  by  M,  and  set  YM 
to  this  quantity.  Set  YMS  to  YMS  + YM,  set  SUM  to  YMS,  and  branch  to  line 
n. 


Line  number 


n + 1 


n + 2 
n + 3 


n+5 


M : N - 1 


M — M + 1 
YM  — YM  * Y/M 


YMS  — YMS  + Y 
SUM  — YMS 
XKS  — 1.  0 
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It  follows  that  if  both  Eq.  (3.  8)  and  IM  + K > M (Y)  are  satisfied  then 


N-l 

oPN<x.Y)=]r 

m^O 


Y Y 


m 


m NIK 


m-N  -X„  3L  k 


_ + Y*  o‘Y  1-  1 1 - V e N -2. 
! m*.  \ Lm  k! 

m=N  \ k=C 


_ m-N  -T9  -rr  k 

Ee-Y  Ym  h V*  e‘XN  XN 

m!  I "La  k! 

=N+K+1  \ k=0 


m=N+K+l 


N-l  N+K  m 

> Y •**  ^ - E e'Y  ^7  Cl  - </*) 

- (—»  m!  a— < m! 

m=0  m=N 


+ Y e-Yl!l 

* m'. 

m=N+K+l 


N+K 

> y e-Y  £1 

- L-t  m! 

N+K 

-i  E 

o 

ii 

£ 

m=N 

> 1 - e/2  - e/2 

= 1 -e 

(3.  12) 


so  to  within  an  accuracy  e , is  equal  to  one.  Second,  problems  of  over- 

flow and  underflow  still  exist.  Underflow  can  be  treated  by  determining  if  X 
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(or  Y)  is  larger  than  a quanta  . E where  e is  at  (or  near)  the  smallest  number 

representable  by  the  computer.  If  X is  greater  than  E,  so  that  evaluating 
-X 

e would  cause  an  underflow,  we  rewrite  the  term 


-X  xk 
e — - e 
k'. 


- X - kf  n X + 


k I 

2>- 

n=l 


= e_Exp 


(3.  13) 


* -Exo 

and  compare  Exp  to  E.  If  Exp  is  greater  than  E,  then  certainly  e H <<  e 

and  can  be  ignored.  We  then  increase  k until  Exp  < E and  start  our  sum- 
mation from  that  value  of  the  index.  The  Y terms  can  be  similarly  handled. 
If  Y is  small  enough  relative  to  N,  the  sum 


-Y  Y 


by  itself  may  be  within  e of  1,  and  this  is  tested  in  the  program. 

If  Y is  large  relative  to  X,  then  QPj^(X,  Y)  is  small  and,  depending 
on  N,  may  be  within  € of  zero.  A test  is  made  in  the  program  (see  Fig.  4, 
KD:X)  to  determine  if  setting  qP^,  to  zero  is  appropriate. 

We  have  assumed  that  the  accuracy  of  the  exponential  function  is 
greater  than  e and  its  errors  can  be  ignored.  It  is  necessary,  however,  to 
consider  round-off  errors.  Since  we  arc  working  in  floating  point,  relative 
errors  are  small  for  a multiplication:  since  all  our  answers  lie  between 

zero  and  one,  the  absolute  error  is  less  than  relative  error  and  can  be  ignored. 
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For  addition,  round-off  errors  can  be  important.  If  the  algorithm  of  Fig.  2 
terminates,  then  the  answer  is  accurate  to  within  e.  Round-off  error,  how- 
ever, can  prevent  it  from  terminating.  The  problem  is  exemplified  by  the 
following  situation:  suppose  XM  is  zero. 


•i-i- 1 „ 

Ee-Y  lH 

6 rn! 


< 1 - e 


so  that  the  criterion,  Eq.  (3.  6),  is  not  met,  and 


where 


-Y  Y , 
e — < 5 


e'Y  IZ]  10‘ 

m!  I 


and  n is  the  number  of  significant  digits.  Then,  because  of  round-off  error. 


L 

V -Y  Yr 
> e — 


E-Y  Ym 
6 ml 

as  represented  m=0 
by  the  computer 


I as  represented 
| by  the  computer 
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so  that  the  algorithm  will  not  end  since  the  criterion,  Eq.  (3.  6),  will  never 
be  met.  The  larger  the  value  of  Y,  the  smaller  each  individual  term  in  the 
sum  and  the  greater  the  number  of  significant  terms.  For  any  e and  any 
finite  number  of  significant  digits,  there  exists  a large  enough  Y such  that 
our  criterion  will  never  be  met. 

We  therefore  have  a trade-off  between  the  number  of  significant  digits 
carried,  the  accuracy,  and  the  maximum  values  for  X and  Y.  If  Z = max(3T,  Y), 
the  worst  case  situation  is  N = 1 and  X = Y = Z.  We  consider  some  specific 
cases.  Withe  - 10  ^ and  double  precision  (16.  8 digits,  or  more  precisely 

56  bits),  the  upper  limit  on  Z was  not  found:  it  is  greater  than  1 million. 

-12 

With  e s 10  and  double  precision,  an  upper  limit  between  110,  000  and 
115,000  exists  on  Z.  For  single  precision  (7.  2 significant  digits  or  24  bits). 
Fig.*  3 shows  approximate  upper  limits  on  Z for  the  different  e.  If  one  can 
be  satisfied  with  the  limits  single  precision  imposes  on  Z for  a given  e,  then 
computational  efficiency  will  be  derived  from  using  it;  otherwise  double 
precision  must  be  used. 

The  final  program  is  shown  in  Fig.  4,  and  the  double  precision 
FORTRAN  version  with  e = 10  and  e = 10  is  given  in  Appendix  A. 

Two  final  points  that  should  be  mentioned  are  (1)  since  it  will  be  use- 
ful in  computing  some  of  the  latter  Pp(s ),  the  value  of  e Y */(N  - 1 ) ! 
is  also  an  output  of  our  program  for  qFj^(X,  Y).  For  values  of  Y for  which 
gPj^fX,  Y)  is  set  to  1 or  0,  we  choose  e~Y  Y^"  */(N-l)!  = 0,  and  (2)  the 
relationship  between  ^P^X,  Y)  and  the  Q - function,  which  is  def.  xed  by 


}Nfa.P)=  f vg)N'* 

Jb 


2,2 

a + v 


Ij^jforv)  dv 


(3.  14) 


Fig.  7>.  Approximate  limits  on  X and  Y vse,  for  single  precision. 


j 


pr? 


we  need  not  calculate  the  second  term  for  N S 2.  The  program  for  ^(X,  Y) 

is  shown  in  Fig.  5,  and  FORTRAN  versions  for  e = 10"^  and  10~^  are  given 
in  Appendix  A. 


3.  3 Case  2: 

In  terms  of  qP^(X,  Y),  2^N^'  Y)  becomes  simply 


JN,x*Y,s»PN(°T7f) 


(3.18) 


3.  4 Case  3: 

Rewriting  Eq.  (2.  6)  in  terms  of  0PN(X,  Y)  we  obtain 


1 + X/2  / , , (X/2)  Y \ 

\ (1  + X/2)2/ 


for  N = 1 


3PN(X,  Y)  =< 


yn  - ie  -y  x 


<N-2>;  (T+Xn/2)  +0PN-  1<°’Y>  for  N > 2 


*('  v) 


N - 2 1 + Xn/2 


I-iLll  +_y 


Xn/2  14XN/2j 


1 ‘ °Pn  - 1 (°' rrtrr)] 


(3.  19) 
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since  Oil-  nP.T  ( 0, = \ 51,  v/e  first  check  to  see  if 

0 N V i+ V2/ 


— 
[l  + xy  2 


- (N  - 


2)  In  4 +^i— -j| 

\ XN/2yj  r ^ N . 2 + Y I 

L'V2  1+V2J 


(3.  20) 


in  which  case  the  last  term  for  N £ 2 can  be  ignored.  As  mentioned  above, 

we  save  the  term  e Y Y'^1  2*  * in  calculating  qP^  ■ 1^’ Y^’  and 

this  is  used  for  the  first  term.  The  program  is  shown  in  Fig.  6,  and  FORTRAN 

*•  6>  1 2 

ver  lions  for  e = 10  and  10  are  given  in  Appendix  A. 


3.  5 Case  4: 


Equation  (2. 7)  i>  manipulated  by  interchanging  the  order  oi  summation 


nRN-I 

-HH  V 

ll  + */*/  li=i< 


1 + X/2  G +X/2)  N!(f) 

ml  “ k!  (N  - k)  ! 


E • 1 


X / Y \m  r—v  k 1 

+ x/2  (1  + x/z)  y»  N{|) 

m!  Z-i  k!(N-k)! 

k— n 


m=2N 


m*.  <k!(N-k)! 

k=0 


M=-f  [y  e ' 1 +Wz  (1  /x/2)  y>  »'(r)k 

V + X/y  mi  k:  (N  - k)  ! 


+/I+*fv 


Y / Y \m"[ 
1+X/2U+X/2  / 


1 +X/2  Vl  +X 


1 +X/2  \1  +X/2 


i k.  (N-k):y1+xyo/  \i+x/2 


There  are  two  cases  in  which  we  can  simplify-  the  calculation  of  P (3?,  Y) 

The  first  is  the  case  where  Y/(l  + X/2)  is  small  enough,  relative  to  N,  that 
0PN  (°'  lTk/zj  **  within  eo£1>  and  we  set  4Pn(3T,  Y)  to  1 ignoring  the  second 

y 

term.  The  second  case  is  when  — — — is  so  large,  relative  to  N,  that 

l+X/2 


2N-1 


X>1  + 


X/2  \l  + X/2) 


and  we  can  set  Y)  to  zero.  This  translates  into  a test  of  the  form 


2N  - 


KK  /-X.-N 

C \l  + X/2/ 


(3.  22) 


where  ^ ^ same  function  as  in  Eq.  (3.  8)  and  depends  on  e.  The 

program  for  ^P^(X,  Y)  is  shown  in  Fig.  7,  and  FORTRAN  versions  for  e = 10~^ 


and  10  4*  are  given  in  Appendix  A. 
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SECTION  4 

APPROXIMATION  FOR  CASE  5,  LOG-NORMAL  FLUCTUATION 


In  this  section  an  approximation  is  presented  for  ^.P^fX,  Y,p  ).  Kramer 
et  aL  , [Ref.  12]  suggested  in  their  Appendix  B that  for  N = 1,  5^N(X,  Y,  p ) 
can  b.i  approximated  by 


.P.{X,  Y,p  ) = -i-  erfc 


(4.  1) 


erfc  (X)  s 1 


/'X 

-4-  f ** 
Jo 


is  the  complementary  error  function.  No  mention  of  accuracy  is  given,  and 
it  became  apparent  that  it  is  dependent  mainly  on  the  parameters  Y,  X,  and  p , 
the  accuracy  deteriorating  if  Y is  large  or  if  X/p  is  near  Y. 

We  will  generalize  the  result  for  any  N and  consider  the  accuracy 
problem.  One  way  of  looking  at  the  approximation  is  to  consider  that  we  have 


■ * * * •*** 


replaced  QPN(X,  Y)  by  zero  for  X less  than  some  and  by  one  for  X greater 
then  Xj.  Then, 


:>NPC,Y,p)=  / px(x|x,p)  dx^ICXj) 

Jx. 


(4.2) 


If  we  define  Ij(Xj)  and  ^(Xj ) as 


„ X 

:«>=/ 

Jq 


Il(Xl)s|  oPN(X,Y)pX{xlX’P)  ^ 


(4.3) 


V'f 

Jx 


*2^X1^  / t1  - 0PN^X*  ^0  PY(x|x,p)dx 


(4.  4) 


then  our  error,  E_ , is 

D 


e5  ■ w - yx.)  . 


(4.5) 


QPN(Xf  Y)  is  monotonic  increasing  in  X from  aero  to  one  and,  since  we 
are  approximating  this  function  by  zero  for  X <X^  and  by  one  for  X <X^, 


we  want 


opn(xi'y>  * 7 


(4.6) 


NXj  would  then  be  the  median  of  the  distribution  density 


0PN(*|V,=  £ e-Y  £ 


-X  X 


m - N 


m=N 


e 0<x<» 

m!  (m  - N)! 


(4.7) 


but  would  also  be  very  cumbersome  to  determine.  Instead,  we  use  the  mean 
ofEq.  (4,  7)  for  NXj 


OO  /•» 

m=N  J 0 


-x  m - N + 1 
e x 


(m  - N)‘. 


QO 

YN  V ym(m  + I),> 
+—J  (N  + m)!  m! 


00 

E(m  + 1 - 


m=N 


e~ Y Y Y*  y—  — 1 - (N  - 1)  e‘Y  Y* 

Lu  (m  - 1)!  m! 

m=N  m=N 


i 

i 





for  N = 1 


Sr 


V'  * 


Y - (N  - 1)  +e 


-Y 


N-2 


r vm  y<n  ’ X>1 

<H-l,YLlr+^nnJtorNfc2 

L m=0  J 


for  N = 1 


[Y  - (N  - I))  [l  - PpA]  + e'Y 


N 


for  N S 2 


(4.8) 


where  Pp^,  the  probability  of  false  alarm,  is  defined  by  Equation  (5.  1). 
V/hen  P„ , is  very  small,  which  is  usually  the  case,  we  can  ignore  it  as  well 

x A 

as  the  term 


-Y  Y 


N 


(N  - 1)! 


and  obtain  as  our  value  for  X, 


x _ Y - (N  - 1) 


N 


(4.9) 


30 


tea 
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NXj  is  near  the  median  and  it  is  shown  in  Appendix  C that  — ■(— . T— J , Y^ 

approaches  1/2  as  Y approaches  infinity;  i.  e.  , as  Y-*  »,  NX^  becomes  the 
median.  Using  Eq.  (4.  9),  the  generalization  of  Eq.  (4.  1)  becomes 


yx,Y(P)=i|. 


/y  - (N  - 1] 

j = i erfc  1 

Y - (N  - 1)\ 

I " NM  1 

\ N 

1 J er-'-  1 

\ V2<r  / 

(4.  10) 


Since  our  approximation 


0PN(x,y> 


(Ofor 


(4.  11) 


improves  rapidly  with  increasing  N,  then  our  approximation  for  ^P^(X.  Y,  p) 
improves  rapidly  with  increasing  N.  The  worst  case  is,  therefore,  for 
N = 1. 

We  now  consider  the  error  in  this  approximation.  We  have 


5PN(X,Y,p)  = KXj)  +1^)  - MXj) 


(4.  12) 


and,  with  Xj  chosen  as  in  Eq.  (4.  9),  Ij  and  are  not  only  small  for  most 
parameter  values  but  also  about  equal  and  tend  to  cancel  out  each  other.  For 
most  cases  we  have 


< 0,  005 


(4.  13) 


SECTION  5 


THRESHOLD  LEVEL 


The  threshold  level,  Y,  is  related  to  the  probability  of  false  alarm, 


PFA>  hV 


PFA  = e 


E If  =0PN<°-y» 


where  N is  the  number  of  incoherent  integrations.  DeLong  and  Hofstetter 
[Ref.  13]  suggest  solving  for  Y in  Ea.  (5.  1)  for  a given  by  solving  for 
the  root  of  f(Y)  where 


f(Y)  = In 


opN(°.Y) 


Employing  the  Newton-Raphson  technique  [Ref.  14]  and  noting  that 


f'(Y)  = - 


rN  - 1 


ON1’ 


(5.  3) 


S 
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SECTION  6 


GENERALIZED  FORMULATION 


The  results  of  Cases  0 to  4 have  been  generalized  and  the  notation 
unified  by  Swerling  [Refs.  15,  16]  by  the  introduction  of  parameter  K.  If 
the  chi-square  distribution  is  generalized  to  a gamma  distribution  with 
parameters  K and  X 


(?|k,X)  = - /ii)e‘K^X  0 < K < 

* 1 r(K)  lx/ 


0 <K  <1 


K = 1 
K = 2 


K = N 


K = 2N 


corresponds  to 

corresponds  to 
c-'-'-responds  to 

corresponds  to 

corresponds  to 

corresponds  to 


Weinstock  Case  [Ref.  9] 
Case  1 


Case  3 


Case  2 


Case  4 


The  probability  of  detection  [Ref.  16]  becomes 


/ K ' 

lK/  XN 

\b 

\k  + xN; 

) \K  + xN 

) 

N-l+b 


-Y  Y 


(6.  1) 


and  Eq.  (6.  1)  can  be  manipulated  as  was  Eq.  (2.  3)  to  obtain 


m-N 


pn,x.  v.k,  Ee_Y  y - e h 


m=N 


i \K  k 


1 K / V'hfJ 


5i  \ 


(6.2) 


A program  to  evaluate  Eq.  (6.  2)  to  within  an  accuracy  e is  given  in  Fig.  8. 

This  program  can  be  modified  so  as  to  eliminate  underflow- overflow  problems 
as  before  by  combining  elements  of  Figs.  4 and  8 so  that  the  variable  K is 
introduced  and  the  quantities 


i— fand/-  ^ 


1 + VK  7 B 


N , ' XN> 


replace  e and  e 


respectively,  in  Fig.  4. 


•M«vtVWUMlV. 
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SECTION  7 
CONCLUSIONS 

Highly  efficient  and  accurate  algorithms  have  been  presented  for 
calculating  the  probability  of  detection  for  Marcum  (Case  0)  and  Swerling 
(Cases  1 through  4)  models.  Tests  and  reformulations  are  included  to  avoid 
the  underflow-overflow  problem  usually  encountered  for  large  parameter 
values. 

A simple  approximation  has  been  provided  for  the  case  in  which  the 
radar  cross  section  is  constant  with  a scan  but  fluctuates  log-normally  from 
scan  to  scan.  The  accuracy  is  limited  but  will  suffice  for  most  applications. 
Guidelines  as  to  which  parameter  combinations  might  pose  an  accuracy  prob- 
lem were  given. 

An  efficient  algorithm  for  determining  the  threshold  level  Y for  a given 

P_  . has  also  been  included. 

FA 

Finally,  with  PpA  =1°  » Figs.  9,  10,  and  11  compared  PD(s)  for 

several  of  the  cases  for  N = 1,  10,  and  100, respectively.  Case  4 curves, 
omitted  on  Figs.  10  and  11,  lie  between  Case  2 and  Case  0.  For  Case  5, 
p = 1.  5 is  used  throughout. 
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PROBABILITY  or  DETECTION 


!0.  Probability  of  detection  vs  for  N 
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APPENDIX  A 


FORTRAN  PROGRAMS 


In  this  Appendix,  FORTRAN  versions  of  the  programs  to  evaluate  the 
probabilities  for  cases  0-4  for  e = 10  and  e = 10  are  given.  In  addition, 
a program  to  determine  the  threshold  for  a given  Pp  . is  given. 

The  subroutine  program  for  qP^{X,  Y)  i3  titled,  with  e = 10-^, 
PNXYS(N,  X,  Y,  P,  YMO)  and,  with  € = 10“12,  PNXYT(N,  X,  Y,  P,  YMO)  where 
N = number  of  incoherently  integrated  pulses 
= average  single-pulse  3ignal -to -noise  ratio 
- threshold  level 


X 

Y 

P 


= 0PN<X'Y> 


YMO  = e 


= «-Y  r 


N - 1 


(N  - 1)! 


N,  X,  Y,  are  inputs;  P and  YMO  are  outputs.  The  output  YMO  is  included 

since  other  programs  calling  on  Qp^fX,  Y)  have  use  for  it. 

Subroutine  programs  for  Swerling  (Cases  l-4,  are  given  in  a single 

program  with  four  entry  points.  The  titles  for  the  four  cases,  with  10  are 

SSWC1(N,  XBAR.  Y,  P),  SSWC2(N,  XBAR,  Y,  P),  SSWC3(N,  XBAR,  Y,  P)  and 

-12 

SSWC4(N,  XBAR,  Y , P),  respectively,  and  those  for  e = 10  are  entitled 
TSWC 1 (N,  XBAR,  Y,  P),  TSWC2(N,  XBAR,  Y,  P),  TSWC3(N,  XBAR,  Y,  P)  and 
TSW4(N,  XBAR,  Y,  P).  In  all  cases  we  have 
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N = number  of  incoherently  integrated  pulses 

XBAR  = average  singie-pulse  signai-to-noise  ratio 

Y = threshold  level 

P = probability  of  detection. 


A subroutine  entitled  PREAM  (X,  CE,  XI,  X2,  X3,  X4,  XK,  IND)  is 
included  as  it  is  called  by  ^P^fX,  Y)  and  ^P^  (X,  Y).  It  evaluates  the  functions 
Ke(X)  and  M£(Y)  and 

X = the  value  of  X or  Y 

CE  = the  corresponding  of  (B:  1 1 ) 

XI  = 18  or  34  for  Kq(X),  and  40  or  80  for  M^(Y) 

X2  = 150  or  175  for  KJX) 

X3  = 0.  75  or  0.  8 for  K (X),  and  1.  3 or  1.  28  for  M (Y) 

£ £ 

X4  = 20  or  55  for  K (X),  and  24  or  52  for  iVL(Y) 

€ c 

IND  is  a flag,  if  IND  = 2,  the  routine  calculator  M^(Y), 
otherwise  K^(X) 

XK  = output,  either  K^(X)  or  M^(Y) 

Finally,  the  subroutine  program  THRESH  (N,  TOL,  PFA,  Yl)  is  used 
to  determine  the  threshold  Yl  in  accordance  with  Eqs.  (5.  4)  and  (5.  5)  where 


N = number  of  incoherently  integrated  pulses 

TOL  = tolerance,  i.  e.  , acceptable  difference  between  , , and 

X\  T 1 

Y^  for  termination  of  the  iterative  process 
Pp^  = input  false  alarm  probability 
Yl  = output  threshold  value. 


A table  of  values  (Table  A-I)  is  included  for  the  purpose  of  checking  out 


the  programs. 


USERID  MASR 
OREAD  PNXYS 


CLASS  A NAME  PNXYS  FORTRAN  02/14/75 

FORTRAN  A1  NASR  8/20/74  14010 

SUBROUTINE  PNXYS(N,X*Y*P( YMO» 

IMPLICIT  REAL*8( A-HtO-2) 

DATA  EtEPSILN/175.D0»l.D-6/ 

X*DFLOATIN)*X 

CALL  PREAMBIX, 16. 1200,18. DO, 150. DO,. 75DO,20.DO,EKX,1> 

CALL  PREAMBI Y, 16. 1200,40.00, 40.00, 1.3D0,24.00,EKY,2» 
IF(DFLOAT(N).LE.EKY-EXX)GO  TO  300 
P*1.00 
YMO*O.DO 
RETURN 
300  I8*N-1 

K*0 

IF(X.GT.E)GO  TO  10 
XX*DEXPC-X1 
5 XKS*XK 
M*0 

IFCY.GT.E1G0  TO  20 
YN*DEXPC-YI 
70  YNS*YM 

40  IFOt.EQ.IBIGO  To  30  . 

M«M*1 

YN*YM*Y/FLOAT ( M ) 

VM’s-f  YM 

IF( EPSILN.LE. l.DO-YMS)GO  TO  40 

P*YMS 

YMO*YM 

X=*X/OFLOAT(N) 

RETURN 
30  SUN=YHS 
YMO=YM 
60  K*Mi 
M=H*1 

YM*YH*Y/FLOAT (Ml 
YMS=YMS*YM 

91  SUM*SUM*YH*{ l.DO-XKS ) 

XX*XK*X/FLQAT IK) 

XKS=XKS*XK 

ANS*( 1.00-YNS)*{ l.DO-XKS ) 

IF( ANS.G1 .EPS II N) GO  TO  60 
50  P*SUH 

X*X/0FL0AT(N) 

RETURN 

20  YHLY=0.00 

YLOG*OLOG( Y ) 

600  M*M+1 

XMF=DLOG( DFLOATCH J ) 

YMLY=YMLY*YLOG-XHF 

IFl Y-YMLY.GT.EJGO  TO  600 

YM*DEXP ( -Y+YMLY ) 

IF(N.L6. I8)G0  TO  70 
YN0=0. 

YMS=YM 

SUM=0.00 
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12028023 


KD=M-N 

IF{ DFLOAT (KD).L7.X)GO  TO  80 
P=0.00 

X=X/DFLOAT(N) 

RETURN 
80  K=K+ 1 

XK=XK*X/FL0A7 ( K ) 

XKS=XKS+XK 
IF(K-KD)80, 90,90 
90  K=K4l 
GO  TO  91 
10  XKLX=0.00 

XLOG=DLOG( X ) 

210  K=K+1 

XKF=OLOG(OFLOAT(K) ) 
XKLX=XKLX+XLOG-XKF 
IF(X-XKLX.GT,E)GO  TO  210 
XK=DEXP{-X4XKLX) 

IB=IB*K 
GO  TO  5 
END 
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USERID  MASR  CLASS  A HAKE  PNXYT  FORTRAN  02/14/75 

OREAO  PNXYT  FORTRAN  A1  NASR  8/20/74  14011 

SUBROUTINE  PNXYTIN,X,Y,P,YM01 
IMPLICIT  REAL*8IA-HP0-2) 

OATA  E.EPSILN/175, DO, 1.0-12/ 

X*DFLOAT (N)*X 

CALL  PREAMBIX, 30.00,34. DO, 175. D0». 800, 55.00, EKX,1) 

CALL  PREAHBIY, 30.00,80.00,80.00, 1.28D0,52.D0»EKY*2) 
IF(DFLOAT(N).LE.£KY-EKX)GQ  TO  300 
P=1.D0 
YN0=0.D0 
RETURN 
300  IB=N-i 
K=0 

IF(X.GT.E)GO  TO  10 
XK=OEXP(~X) 

5 XKS=XK 
H«0 

IF(Y.GT.E)GO  TO  20 
YN=DEXP(— Y) 

70  YHS=YH 

40  IF(M.EO.IB)GO  TO  30 
N=H*1 

YN*YN*Y/FLOAT(M) 

YNS-YHS+YM 

IF( EPSILN.LE. l.DO-YMS )G0  TO  40 

P=YMS 

YMO=YN 

X=X/DFLOAT ( N) 

RETURN 
30  SUM=YMS 
YMO=YM 
60  K=K+1 

YH=YM*Y/FLOAT(H> 

YNS=YHS*YM 

71  SUM=SUM*YM*( l.DO-XKS) 

XK=XK*X/FLOAT ( K } 

XKS=XKS*XK 

ANS=( l.D0-YHS)*( l.DO-XKS) 

IF CANS.GT.EPS I LN) GO  TO  60 
50  P=SUN 

X=X/DFLOAT (N) 

RETURN 

20  YNLY=0.D0 

YLOG=OLOG( Y ) 

600  N=N+1 

XNF=DLOG( OFLOAT(M) ) 

YNLY=YMLYvYLOG-XHF 
!F(Y-YHLY.GT.E)GO  TO  600 
YM=DEXP ( -Y+ YML Y ) 

IFIH.LE. IB )G0  TO  70 
YH0=0. 

YMS=YM 

SUN=0.D0 
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I 


KD=M-N  | 

IF  COFLOAT (KD).LT.X)GQ  TO  80  I 

P=O.DO  I 

X*X/DFLOAT ( N)  I 

RETURN  1 

80  K=K*1  | 

XK=XX*X /FLOAT  IK)  S 

XKS*XKS*XK  I 

IF(K~K0)80,90,90  1 

90  K=X*1  I 

GO  TO  91  J 

10  XKLX=0.00  I 

XLOG=CU.OG(X)  i 

210  K=K4-1  a 

XKF=OLOG(  OFLOAT { K ) ) 5 

XKLX=XKLX*XLOG-XKF  I 

IF(X-XKLX.GT.E)GO  TO  210  I 

XK=DEXP  ( -X+-XKLX ) I 

I B= I B+K  | 

GO  TO  5 I 

ENO  I 

> 

1 


I 

¥ 

1 

$ 

| 


I 


■g- 
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USERID  NASR  CLASS  A NAME  SSCASE  FORTRAN  02/14/75 

OREAO  SSCASE  FORTRAN  A1  NASR  8/22/74  12049 

SUBROUTINE  SSWC1IN, XBAR,Y,P> 

IMPLICIT  REAL  *8U-H,  0-2) 

DATA  E, EPS I LN/ 175. 00, l.D-6/ 

XBAR=OFLOAT ( N )*XBAR 
CONST=Y/( 1.D0+X8AR) 

IFIN.GT o 1)G0  TO  5 
IF(CONST.GT.E)GO  TO  6 
P=DEXP( -CONST ) 

GO  TO  100 
6 P=0.D0 
GO  TO  100 
5 SUM=0.D0 

EX=C0NST-DFL0AT(N-1 ) *DLOG(  l.DO+l.DO/XBAR) 

IF( EX. GE. 30.00)50  TO  10 

CALL  PNXYS(N-1,0.D0,Y/I l.+l./XBAR) ,??, YSUH) 

SUM=( l.DO-PP)*OEXP(-EX) 

10  CALL  PNXYS(N-1»0.D0»Y»P1»YSUN) 

P=SUH+P1 

100  XBAR=XBAR/DFLOAT  f N ) 

RETURN 

ENTRY  SSKC2IN,XBAR,Y,P) 

CONST=Y/(  l.DO+XBAR) 

CALL  PNXYSI N» 0.00, CONST, P» YSUH) 

RETURN 

ENTRY  SSWC3IN,XBAR,Y,P) 

XBAR=DFLOAT IN)*XBAR 

C3=1.00+X6AP./2.00 

C0NST=Y/C3 

CZ* 1.D0+2. DO/XBAR 

IF  < N. EQ. 1 )G0  TO  20 

SUK=0.D0 

EX-CONST-OFLOAT ( N-2 ) *DLOGC  C2 ) 

C=Cf€XP(-EXI *(  l.D0-2.D0*DFL0AT(N-2) /X0AR+Y/C3) 

IF; C.LE.EPS ILN) GO  TO  200 

CALL  PNXYSI N-l, O.DO»Y/C2» PI » YSUH) 

SUM=C*( 1. DO-Pi) 

200  CALL  PNXYSIN-1,0.D0,Y,P2» YSUH) 

SUM=SUM+P2 
P=SUM+YSUM* 1. D0/C3*Y 
GO  TO  250 

20  IFICONST.LT. E)GO  TO  21 
P-O.DO 

GO  TO  250 

21  P=DEXPI -CONST ) * l l.DO+XBAR/ ( 2. D0*C3 ) *Y/C3 ) 

250  XBAR=XBAR/DFLOAT (N) 

RETURN 

ENTRY  SSHC4 1 N * XBAR  , Y , P ) 

C3= l.DO+XBAR/ 2. DO 
XN*0FL0AT  CN)*OLOG(C3) 

CONST  =Y/C3 

CALL  PNXYSI N.O.DO, CONST, SUM, YSUH) 

I FI l.DO-SUH.GT.EPSILN)GO  TO  30 
P=1.00 


1202902A 


50 


GO  TO  300 

30  CALL  PREAHB ( CONST t 16* 12D0, 18. DO » 150. DO, . 75D0 ,20. DO » 6KY» 1 ) 

IP( OFLOAT ( 2*N-1 ) .GE.EKY1G0  TO  31 

P=0.00 
GO  TO  300 

31  M=N-1 

IF(XN.GT.E)GO  (0  40 
2K»DEXP(-XN) 

51  2XS=ZK 

Z*XBAR/2.D0 

42  IF(N.EQ.2*N-l IGO  TO  41 
M*M*1 

YSUH=YSUM*CONST/OFLOAT(H) 

SUM=SUM*YSUN*< 1.00-ZKSI 

ZK=ZK*Z*( 2«D0*DFL0AT ( N ) -DFLOATtN J ) / (DFLOAT (M)-DFLOAT ( N) ♦ l.DQ) 
ZKS=ZKS*ZK 
GO  TO  42 
41  P=SUH 
300  RETURN 
40  H=M-»1 

XN=XN-DLOG(XBAR/2.DO»-DLOGC2.DO*DFLCAT(N)-DFLOAT(M) )/<DFLOAT(M)-DF 
1L0AT (N)  + 1 ) 

IFtXN.GT.E )G0  TO  50 
ZK=DEXP(-XN) 

GO  TO  51 

50  YSUH=YSUN*CONST/OFLOAT (M) 

SUM=SUM+YSUM 
GO  TO  40 
END 
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USERID  NASR  CLASS  A NAME  TSCASE  FORTRAN  02/1A/75 

OREAD  TSCASE  FORTRAN  A1  MASR  8/29/74  11014 

SUBROUTINE  TSWCHN.XBAR, Y,P) 

IMPLICIT  REAL*8( A-H,0-Z ) 

DATA  E, EPS ILN/ 175.00, i.D-1 2/ 

X8AK= DFLOAT ( N 1 *XBAR 
CONST* Y/( l.DO+XBAR ) 

IF(N.GT.1)G0  TO  5 
IF (CONST .GT.E )G0  TO  6 
P*DEXP(-CONST) 

GO  TO  100 
6 P=O.DO 
GO  TO  100 
5 SUM=O.DO 

EX*CONST-DFLOAT (N-l  )*DLOG(  l.DO+l.DO/XBAR) 

IF{ EX. GE. 30.00 )G0  TO  10 

CALL  PNXYT ( N-l, 0.00, Y/( l.+l./XBAR) , PP»YSUM) 
SUM=(1.00-PP)*0EXP(-£X) 

10  CALL  PNXYT(N-1,0.D0,Y,P1,YSUM) 

P=SUM*P1 

100  XBAR=XBAR/DFLOAT(N) 

RETURN 

ENTRY  TSWC2(N,XBAR,Y,P) 

C0NST=Y/(1.D0*XBAR) 

CALL  PNXYT ( N, 0. DO, CONST » P» YSUM 1 
RETURN 

ENTRY  TSWC3(N,XBAR,Y,P) 

XBAR=DFLOAT(N)*XBAR 

C3=1.D0*XBAR/2.D0 

C0NST=Y/C3 

C2-1.D0+2.D0/XBAR 

IFlN.EQ.l JGO  TO  20 

SUM*C.00 

EX=CONST-OFLOAT(N-2)*OLOG(C2> 

C=OEXP(-EXi *1 l.00-2.D0*DFL0AT(N-2)/XBAR*Y/C3) 

IF (C.LE.EPS ILN) GO  TO  200 

CALL  PNXYT l N-l, 0. DO, Y/C2, PI, YSUM) 

SUM=C»( i.OO-Pi) 

200  CALL  PNXYT ( N-l, 0. DO, Y, P2, YSUM) 

SUM=SUM»P2 

P*SUM+YSUM^1.00/C3*Y 
GO  TO  250 

20  IF (CONST.LT .E)GO  TO  21 
P=0.D0 

GO  TO  250 

21  P=DEXP( -CONST )♦( 1.00+XBAR/ ( 2.D0*C3 ) *Y/C3 ) 

250  X8AR=XBAR/DFL0AT( N ) 

RETURN 

ENTRY  TSWC4(N,XBAR,Y,P) 

C3=i.00*X8AR/2.D0 

XN--0FL0AT(N)*DL0G(C3) 

C0NST=Y/C3 

CALL  PNXYT(N,0. DO, CONST, SUM, YSUM) 

IF ( 1. DO-SUM. GT.EPSILNJGO  TO  30 
P=1.D0 


Si?' 


E£ 


GO  TO  300 

30  CALL  PREANBICONST, 30.00,34.00, 175.D0,. 800. 55.00, EKY.l) 

IF! DFLOAT ( 2*N- 1 ) .GE • EK Y ) GO  TO  31 

P*0.00 
GO  TO  300 

31  M*N-1 

IFlXN.GT.EIGO  TO  40 
ZK*DEXP(-XNI 
51  ZKS=ZK 

l-X BAR/2.00 

42  IF!M.EQ.2*N-1)G0  TO  41 
H=H+i 

YSUM=YSUM*CONST /DFLOAT t M ) 

SUM*SUM+YSUN*! l.DO-ZKS) 

ZK*ZK*Z*( 2. DO*DFLOAT ( N I-DFLOATl NJ ) / (DFLOAT !M)-DFL0AT!N)+1.D0) 

IF! MODI M»9) .NE.O )G0  TO  47 
WRITE!6,48)ZK 
48  FORMAT! IX, ‘ZK**, 015. 8) 

47  ZKS*ZKS*ZK 
GO  TO  42 
41  P=SUM 
300  RETURN 
40  M*M*1 

XN*XN_OLOG! XBAR/2.DO)~OLOG!2.DO*DFLOAT! NJ-DFL0AT1M) l/IDFL0AT(M)-0F 
1L0AT(N)+1 ) 

IFlXN.GT.EIGO  TO  50 
ZK=OEXP(-XN I 
GO  TO  51 

50  YSUM=YSUM*CONST /DFLOAT ! M > 

SUM=SUM*-YSUM 
GO  TO  40 
END 
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USERID  HASR  CLASS  A HAKE  PREANB  FORTRAN 

OREAD  PREANB  FORTRAN  A1  HASR  8/20/74  150X6 

SUBROUT  I NE  PRE AMB( X « CE  » X 1 * X2 , X3 , X* , XK, I ND1 
IMPLICIT  REAL*8(A-H,0-2) 

DATA  TWOP 1/6. 283 18 5 300/ 

IF ( IN0.EQ.2iG0  TO  20 
IFU.GT.XUGO  TO  1 
XK*0.00 
RETURN 

1 IF(X.GT.X1.AND.X.LT.X2)G0  TO  21 

XK*X3*X-X4 
RETURN 

21  CONST— 1.00 
2 XK*X*CONST*DSQRT  C CE*X ) 

5 FX=X-XK+XK*0L0G( XK/X 5*DL0G( TWOPI ♦XK l-CE 
FXP=DLOG{ XK/X  J + .5D0/XK 
XK 1=KK-FX/FKP 

I F ( DABS ( XK 1-XK ) .LT. . IDO ) RE  TURN 
XK-XK1 
GO  TO  5 
20  C0NST=1.D0 

I F ( X.LT.X1 ) GO  TO  2 

XK*X3*X+X4 

RETURN 

END 


02/14/75  12028008 


S 5 
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USERID  MASR  CLASS  A NAME  THRES  FORTRAN  02/14/75  12028014 

OREAD  THRES  FORTRAN  Ai  MASR  8/22/74  15025 

SUBROUTINE  THRESH!N,T0L,PFA,Y1) 

IMPLICIT  REAL*8I A-HiO-Z) 

OATA  E/175.00/ 

PL0G*DL0G10!PFA) 

XPF A* 1 • OO/PFA 
IF1N.GE .40.00) GO  TO  30 
Y=DFL0AT1N) **1. 15-2.DO*PLOG 
GO  TO  31 

30  Y=DFLOAT! N)-8.DO*PLOG 

31  IFIY.GT.EIGO  TO  10 
212  M*0 

YLX=0.D0 
YLOG-OLOG! Y) 

TSUM=DEXP(-Y) 

21  IF1M.EQ.N-DG0  TO  20 
110  M=M*1 

SUHK=DLOGI DFLOAT(M) ) 

YLX=YLX*YLOG-SUMK 

TSUM=TSUM*DEXP!-Y*YLX) 

GO  TO  21 

20  Y1=Y*DL0G f XPF A*TSUM ) *TSUM/QEXP ( -Y*YLX I 
IFCDABSI Y1-Y).LE.T0L)RETURN 
Y=Y1 

GO  TO  31 

10  M=0 
YLX=0.D0 
YLOG=DLOG( Y ) 

11  M*M*i 

SUMK=OLOG{OFLOATCM) ) 

YLX=YLX*YLOG-SUMK 
IF( Y-YLX.GT . E )G0  TO  11 
TSUM=DEXP  f -Y+YLX ) 

IF! M-N+l ) 1 10* 20 1 20 
END 

IMPLICIT  REAL*8!A-G,0-Z) 

1000  READ! 5t*)N«T0LtPFA* IEND 
CALL  THRESHCN,TOL,PFA,Y) 

NRITEI6* 100 )Y 
100  FORMAT! lXt*Y=,»015.8) 

IF! I END. EQ.O) RETURN 

GO  TO  1000 

END 


APPENDIX  B 

DETERMINATION  OF  THE  FUNCTIONS  K (X)  AND  M (Y) 

€ € 


The  function,  K.(X),  is  defined  to  have  the  property  that  for  all  integers 


K such  that 


K < K (X) 


(B.  1) 


— < c/2 


(B.2) 


K^(X)  = 0 denotes  that  the  sum  is  empty. 

M (Y)  is  defined  to  have  the  property  that  for  ail  integers  N such 


N 2 Mc(Y) 


(B.  3) 
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I IPM  H > *nr  WIJC  W.V  »*N.**’  **  f 


then 


N-l 

E 

m=0 


-Y  Y 


m 


m ! 


> 1 - c/2 


(B.  4 ) 


or  equivalently 


E 

m=N 


-Y  Ym  /? 
e — - < c/2  . 

m! 


(B.  5 ) 


Using  the  square  bracket  notation  to  mean  the  largest  integer  less  than  or 
equal  to  the  quantity  in  the  bracket,  then  ideally,  we  should  determine  K^(X) 

such  that  [K  (X)]  is  the  minimum  integer  for  which  Eq.  (B.  2)  is  true,  and 

€ 

[M  (Y)]  is  the  maximum  integer  for  which  Eq.  (B.  4)  is  true;  it  is  of  little 
practical  significance  to  have  the  exact  values.  Values  close  to  their  mini- 
mum and  maximum,,  respectively,  are  much  easier  to  derive. 

In  order  to  derive  a useful  K^fX),  we  begin  by  noting  that 


< I 


K 


k=0 


k! 


X^  L + K + K(K  - 1)  + + 

K|  L X yZ  jfk 


s e 


-X  X 


K 


K! 


1 + 


MS'*  *©1 


= e 


xf_ 

K! 


i . a 

\X> 


K+r 


1 -*S 


< e 


£ 


WWj 


1 
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K 

X"X 


(B.  6) 


for  K,  where 


for  e = 10  ^ 

I 

-12 

for  e = 10  16 

K 

the  solution  can  be  used  for  K (X)  provided  — is  less  than  4/5.  The  Newton 

^ X 

Raphson  iteration  technique  to  solve  Eq.  (B.  1 1)  is 


116.  12 
30.  0 


K 


n + 1 


fx(Kn.£) 

fX<Kn  •«> 


(b.  i; 


where 

fX<Kn-)=ln@+2^  ' <813 

The  initial  value  for  K is  obtained  for  Eq.  (B.  11)  by  approximating  InZ  by 
Z - 1 and  ignoring  the  term  j £ n2tiK.  Eq.  (B.  11)  then  becomes 

X - K + K 


or 


K2  - 2XK  + X2  - C X * C:  . 


From  this  table  we  have,  with  g = 10  ^ , that  for  X £ 300,  then  — 5 = 

X 300 

0.727  < 0.8;  and  with  g = 10-12,  that  for  X < 700,  then  ~ < = 0.741  < 0.8. 

For  these  values  of  X we  can  use 


i 1 

i SI 


a 


K (X)  = root  of  fv(K,g) 


(B. 16) 


For  larger  values  of  X,  we  need  another  expression.  By  curve  fitting,  we 

obtain  for  K (X) 

£ 


K (X)  = 0.  75X  - 20 


(B. 17) 


for  X > 150.  Since  Eq.  (B.  17)  is  much  more  simple  than  Eq.  (B.  16),  we  use 

- 12 

it  down  to  X > 150.  For  g = 10  , we  similarly  obtain 


K (X)  = 0.  8 X - 55 


(B. 18) 


_ L — 1 7 

for  X > 175.  For  X < 18  with  c = 10  and  for  X < 34  with  g = 10  , we  set 

K (X)  to  zero.  Summarizing  our  results  we  have 


for  X < 34 


K ,(X)  = { root  of  fY(K,  10’6)  for  18  s X < 150  (B.  19) 

io'b  ’ X 


0.75  X - 20 


for  150  < X 
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for  X < 34 


K _1?(X)=  { root  of  fv(K,  10“14)  for  34  < X < 175 . (B.  20) 


for  175  < X 


M (Y)  is  dealt  with  in  a similar  manner  to  that  of  K (X).  If  N > M (Y), 
e £ € 


then  Y/N  < 1,  we  have 


E-Y  Ym  _ -Y  YN  . , Y 
m!  C N!  N + 1 


+ t + 

(N  + I)  (N  + 2)  *“ 


then  Eqs.  (B.  4)  and  (B.  5)  are  satisfied.  If  Y/N  < 4/5,  then 


or  equivalently 


N + N £n{N/X)  + j Jfn  Zv  K 


e/10 


(B.  2 


which  reduces  as  in  £.q.  (B.  11)  to 


fY(N.  e)  = 0 


(B.  2 


An  empirically  determined  table  of  values  of  M^(Y)  are  given  in  Table  B-IL 


Table  B-H 


&35£0j^3& 


Although  Eq.  (B.  24)  is  the  same  as  Eq.  (B.  11)  except  that  we  have  replaced 
X by  Y and  K by  N,  we  are  not  looking  for  the  seme  root.  fy(N,e)  ha3  two 
roots  and  this  time  we  need  the  larger  root;  whereas  for  f^(K,c)  we  needed 
the  smaller.  This  is  accounted  for  by  our  initial  value  for  N,  which  uses  the 
plus  sign  in  Eq.  (B.  14) 


Nn  = Y + 1C  Y 
0 y € 


(B.  25) 


The  root  can  be  used  for  large  values  of  Y,  but  since  we  have  a more  simple 
expression 


M , (Y)  = 1.  3 Y + 24  for  Y > 40 

10 

and 

M . ?(Y)  = 1.  28  Y + 52  for  Y > 80 

10‘A* 

we  use 

I root  of  fN(Y,  10"6)  for  Y < 40 
1.  3 Y + 24  for  Y £ 40 

and 

!rootoffN(Y,10-12)  for  Y < 80 
1.  28  Y + 52  for  Y s 80 
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APPENDIX  C 


ASYMPTOTIC  VALUE  OF  QPN  ^Y— -— ) , Y^ 


We  show  that 


«_  D /Y  - (N  - 1)  v\  _ 1 

0 N\  N ’ ) 2 


Gradshten  and  Ryzhik  [Ref.  17]  (p.  717,  6.  63.  8)  give  us 


ri  2 (”il> 

I zn  e~  yZ  In  _ 1(2yz)  dz  = eY  - e ^ ' ^S2^ 

JQ  k=-(n-l) 


From  Eq.  (2.  1)  we  have 


1 -opn<x'y)  = 


«‘(V  + MX)W* 


(C.  3) 


1 

2 


1 - e-2'Y  - » - 1)1  £ IkPtY  - ,N  - 

k=-(N- 1) 


+ 2[Y  - (N-  1)]  e'tY‘*N' 


I; 


fY-(N-l) 


ZN  e‘IY“tN’i)l  In(2[Y-(N-1)]2)  dZ 


(C.  5) 


Taking  the  limit  as  Y-*» 


yl 


Y 

Y - (N  - 1) 


1 


(C.  6) 


and 


e-2[Y-(N-l)]  Ik(2[Y  -(N-l)])-*0  £or  |k|<N-  1 


(C.7) 


so  that  Eq.  (C.  1)  is  obtained. 


m imfflH  uaMaaM 


USERID  MASR  CLASS  A NAME  PNLN  FORTRAN  02/14/75  12029019 

OREAD  PNLN  FORTRAN  A1  MASR  8/27/74  14039 

SUBROUTINE  PNLN!N»X»Y»RH0,P5) 

XN=X/RHO 

SIG®  SORT (2.  *AL0G(RH0>) 

ARG=ALOG!!Y-  FLOAT(N-l))/(  FLOAT! N>*XM) )/( 1*414213562  ♦SIG) 

CALL  EREC(  ARC*  ERC ) 

P5*.5  *ERC 

RETURN 

END 

SUBROUTINE  ERECiX.EKC) 

DATA  A1,A2, A3, A4,A5,P, PI/. 225836846, -.252128668, 1.25969513, 
1-1.287822453, .94064607, . 3275911,3. 14159265/ 

SPI=SQRT C PI  I 
IF(X) 10,2,2 
2 C1=0.0 


C2*1.0 

3 ETA-l./I l.-*P*X) 

POLY®! ( ! 1 1 A5*ETA+A4I*ETA*A3)*ETA+A2I*ETA*A1 )*ETA) 
ERC=C1*C2*2.*P0LY*EXP!-X**21 /SPI 
RETURN 
10  Cl*2. 

C2*-l. 


X*-X 
GO  TO  3 
END 


s* 


! 1 

! 

I 

£ 

$ 
& 
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